Microscopic theory of weakly coupled superconducting 
multilayers in an external magnetic field 

Sergey V. Kuplevakhsky 
Department of Physics, Kharkov State University, 

310077 Kharkov, Ukraine 
e-mail : Sergey. V . Kuplevakhsky® uni ver . kharkov . ua 
and Institute of Electrical Engineering, 
SAS, 842 39 Bratislava, Slovak Republic 



Abstract 

We present the first fully microscopic, self-consistent, and self-contained the- 
ory of superconducting weakly coupled periodic multilayers with tunnel bar- 
riers in the presence of externally applied parallel magnetic fields, in the lo- 
cal Ginzburg-Landau regime. We solve a nontrivial mathematical problem 
of a microscopic derivation and exact minimization of the free-energy func- 
tional. In the thin-layer limit that corresponds to the domain of validity 
of the phenomenological Lawrence-Doniach model, our physical results strik- 
ingly contrast with those of our predecessors. In particular, we completely 
revise previous calculations of the lower critical field and refute the concept of 
a triangular Josephson vortex lattice. We show that Josephson vortices pen- 
etrate into all the barriers simultaneously and form peculiar structures that 
we term "vortex planes". We calculate the superheating field of the Meiss- 



ner state and predict hysteresis in the magnetization. In the vortex state, 
the magnetization exhibits distinctive oscillatory behavior and jumps due to 
successive penetration of the vortex planes. We prove that the vortex-plane 
penetration and pinning by the edges of the sample cause the Fraunhofer pat- 
tern of the critical Josephson current. We calculate the critical temperature 
and the upper critical field of infinite (along the layers ) multilayers. For 
finite multilayers, we predict a series of first-order phase transitions to the 
normal state and oscillations of the critical temperature versus the applied 
field. Finally, we discuss some theoretical and experimental implications of 
the obtained results. 



PACS numbers: 74.50. +r, 74.80.Dm 

I. INTRODUCTION 

In this paper, we present the first fully microscopic, self-consistent, and self-contained 
theory of superconducting weakly coupled periodic multilayers (superlattices) of the S/I- 
type (S for a superconductor, I for an insulator or a semiconductor) in the presence of 
externally applied parallel magnetic fields, in the local Ginzburg-LandauS (GL) regime (i. 
e., temperatures are close to the bulk critical temperature T c q, all characteristic dimensions 
of the S-layers are much larger than the coherence length £o-) 

Our approach introduces radically new concepts both into the construction of the theory 
and understanding of the underlying physics. 
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In Section II, we derive a microscopic free-energy functional that describes a smooth 
transition from the single-junction case to the thin-layer limit, when the S-layer thickness a 
is small compared to all other relevant length scales. [By contrast, previous treatment was 
based predominantly on the phenomenological Lawrence-Doniachi (LD) model, applicable 
only in the thin-layer limit.] Our analysis of implications of gauge invariance reveals novel 
facets of the theory that have been completely overlooked up to now. In particular, we estab- 
lish for the first time the existence of fundamental constraint relations coupling the phases 
of the order parameter to the vector potential and equalizing the phase differences at neigh- 
boring barriers. Mathematically, these constraints complement the usual Euler-Lagrange 
equations and make the free energy a minimum. Physically, they state that the average 
intralayer currents are always equal to zero and the local magnetic field has, in general, the 
periodicity of the multilayer. (The latter has recently received strong experimental support 
in polarized neutron reflectivity measurements on artificial Nb/Si multilayers.^) As a result 
of this investigation, we obtain a closed, self-consistent set of mean-field equations. These 
equations have a particularly simple form in the thin-layer limit: Remarkably, the phase 
differences (the same at all the barriers) obey a Ferrell-Prange-typei equation with a single 
length scale, the Josephson penetration depth Aj = (87repjo) ^ 2 (p is the period, jo is the 
critical Josephson current). (This should be contrasted with a mathematically ill-defined 
infinite set of differential equations containing two length scales, proposed without appro- 
priate justification in the literature .§0) Furthermore, due to the absence of screening by the 
intralayer currents in the thin-layer limit, the local magnetic field proves to be independent 
of the coordinate normal to the layers. 

In Section III, we obtain exact analytical solutions to the equations of the thin-layer limit. 
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This limit corresponds to the domain of validity of the phenomenological LD model, which 
allows us to draw a comparison with the results of our predecessors. Unfortunately, for the 
lack of a closed, self-consistent set of mean-field equations minimizing the free energy, some of 
the previous studies have rather spurious character. Thus, the above-mentioned fundamental 
constraints of the theory automatically rule out any possibility of suggested single Josephson 
vortex penetrationll! as well as the occurrence! of a triangular Josephson vortex lattice. 
Another example: in Ref. (0) concerned with the Josephson effect, inability to carry out 
self-consistent calculations misled the authors into an absolutely incorrect conclusion that the 
Fraunhofer pattern of the total critical current occurred in the absence of Josephson vortices. 
Here we undertake a complete revision of these issues. However, our consideration contains 
as a limiting case a particular exact solution for the vortex state obtained in the framework 
of the LD model by TheodorakisS and fully supports phenomenological calculationsQ0 of 
the upper critical field H c2oo in infinite (along the layers) multilayers. 

Among new physical results of Sec. Ill, there are the following. We provide the first com- 
prehensive description of the Meissner state in semi-infinite (along the layers) multilayers and 
show that at the field H s = (ep\j)~ l (the superheating field) the Meissner phase becomes 
unstable with regard to Josephson vortex penetration. We predict simultaneous and coher- 
ent penetration into all the barriers. (This prediction has been confirmed experimentally.^) 
We show that in the absence of screening by the intralayer currents (see above) the "tails" 
of Josephson vortices overlap in the layering direction forming peculiar structures, "vortex 
planes". The lower critical field at which the formation of a single vortex plane becomes 
energetically favorable in an infinite multilayer is found to be i? c ioo = 2(7repAj) _1 . For 
the lower critical field in a finite multilayer with W <C Xj (W is the S-layer length) we 
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obtain H c iw = n/epW, which corresponds to the first minimum of the Fraunhofer pat- 
tern. We prove that the Fraunhofer oscillations occur due to successive penetration of the 
vortex planes and their pinning by the edges of the sample. We show that vortex-plane 
penetration leads also to jumps of the magnetization. (Such features have been already 
observed.@) For a certain field range, we predict a small paramagnetic effect. We calculate 
the critical temperature and the upper critical field of an infinite multilayer. The obtained 
implicit dependence if c2oo (T) exhibits the well-known "3D-2D crossover" and is free from 
the unphysical "low-temperature" divergence^! of the LD model. In addition, we predict 
novel size-effects in finite multilayers: a series of first-order phase transitions to the normal 
state and oscillations of the critical temperature versus the applied field. 

In Section IV, we discuss some theoretical and experimental implications of the obtained 
results. In Appendix A, we write down a few mathematical formulas related to the applica- 
tion of Mathieu functions in Sec. III. 



A. A derivation and exact minimization of the microscopic free-energy functional 



Our starting point is a microscopic second-quantized BCS-type Hamiltonian of the 



II. BASIC EQUATIONS OF THE THEORY 



form@0 



H = |d 3 r^(r) 



1 



V - ieK{r)-ieA ext {r) 



2m 



R 



|i J rf 3 r^(r)^+(r)^ a (r)^_ a (r)^ a (r) + j d 3 n/;+(r)VW p (r)^(r) 
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+ U J d 3 r^(r)Mr) + d 3 rh\v), (1) 

Rb Ft 

R s = U^_ 00 R Sn ,R b = Un^-ooRbnjR = R s U R b , 
R Sn = [—a/2 + np < x < a/2 + np] x [L y \ <y< L y2 ] x (—00 < 2 < +00), 
Rb n — [ a /2 + (n — l)p < x < —a/2 + np] x \L yl <y< L y2 ] x (—00 < z < +00), 

h(r) = V x A(r), H = V x A ext (r) = (0, 0, H), (2) 

V = (V„V„V 2 ).(^,|,|). 

Here h = c = 1, £7^? = hp /2m is the Fermi energy (with fc^ being the Fermi momentum), 
R s and Rb correspond respectively to the superconducting and barrier regions (with a being 
the S-layer thickness, b the barrier thickness and p = a + b the period, the x axis being 
normal to the barrier interfaces), i^ Q (r) is the electron field operator for spin a (a summation 
over repeated spin indices is implied), g < is the BCS coupling constant, Vi mp (r) is the 
nonmagnetic impurity potential, Uq > is the repulsive barrier potential, A ext and A are 
the external (classical) and induced (operator) vector potentials.0 The system is taken to 
be infinite in the direction of the x and z axes, while no restrictions on the linear dimensions 
along the y axis is so far imposed. The external magnetic field H is directed along the z 
axis (see Fig. 1). 

Using field-theoretical methods of Ref. (0), we can derive from ([]]) a microscopic free- 
energy functional of the system Q [A n , A*, A; H], where A n and A are classical variables: 
A n is the pair potential (order parameter) of the nth S-layer, and A = A + A ext is the 
total vector potential, h(r) = V x A(r) being the corresponding local magnetic field. For 
external fields satisfying the quasiclassical condition H <C kp/e^o, in the GL regime 
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£o < a, W = L y2 - L wl , (4) 

where T c0 is the bulk critical temperature, £ — v /27iT c0 is the BCS coherence length 
(v = kp/m), this functional takes on the form 



H 2 (T) T +°° a/ r 

Q[f n , (f)n ,A x ,A y ;H] = ^-^W z Jdy £ j 



a/2+np 



dx 



-a/2+np 



+C 2 (T) £ [[V,/ n (x,y)] 2 + [V,0„(x,y) - 2e^(:r,y)] 2 /„ 2 (x,2/) 



2a£ 



J 2 ,! (a/2 + (n - l)p, y) + f 2 (-a/2 + np, j/) 



-2f n (-a/2 + np, y) /„_i (a/2 + (n - y) cos $„, n _i 
+4e 2 C 2 (T)A 2 (T) Jdx[h(x,y)-Hf 



Lxi 



—a/2+np 

®n, n -i(y) = 0n(-«/2 + np,y) - n _i (a/2 + (n - - 2e ^ dxA x (x, y), 

a/2+(n-l)p 



(5) 



3tt 2 



a = 



7C(3)X (WO 



16Ept 2 (U - Ept 2 ) 
D[t) = r75 exp 



/7 2 



JdttD{t), 
o 

-2b^2m (U - Ept 2 ) 



+ 0O 



X (Co/0 = ^ttttE (2n + 1) (2n + 1 + 6/0 , 



H 2 (T)=4nN(0)A 2 oo (T)r, 



h(x,y) = ^^- d -^l. (6) 
In Eq. (El), we have introduced the reduced modulus < f n < 1 and the phase <fi n of the pair 



potential in the nth S-layer via the relation A n = Aoo f n exp (i<f> n ) , where A DO (T) = y~j^^~ 
is the bulk gap, C( m ) is the Riemann zeta function.0 The rest of notation are as follows: 
W z = L Z 2—L z i is the length of the system in the z direction, D(t) is the tunneling probability 
of an insulating barrier between two successive S-layers [D(l) <C 1], \ (Co/0 is the impurity 
factorlH (I is the electron mean free path), £(T) = £o /^§T^ is the GL coherence length, 
A(T) = 4k(fo fg (Q)Tl " 1/2 is the GL penetration depth, N{0) = mk F /2ti 2 is the one-spin 
density of states at the Fermi level, and H C (T) is the bulk thermodynamic critical field 
near T c0 .El1 The term proportional to a C 1 determines the interlayer Josephson coupling. 
Equation (|6]) is merely the Maxwell equation for the local magnetic field h = (0, 0, h). 

The microscopic free-energy functional (|J) covers all well-known limiting cases. In the 
limit a = (no Josephson interlayer coupling), equation (||) reduces to a sum of free- 
energy functionals of independent S-layers. Making a shift of the coordinate system x — > 
x — a/2 — 6/2 and taking the limit a — >• oo, one gets the case of a single SIS junction. 
Shifting x — ■> x — a/2 and taking a — > oo, b — ► oo, we recover the limit of a semi-infinite 
superconductor in contact with vacuum. 

Our task now is to establish mean-field equations of the theory, which is mathematically 
equivalent to the problem of minimization of @ with respect to /„, <f> n , and A = (A x , A y ,Q). 
This problem should be approached with certain caution, because Euler-Lagrange equations 
for (fi n , and A x , A y are not independent. 
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Indeed, the functional (|^) is invariant under the general gauge transformation 

<f> n (x, y) -> (j) n (x, y) + r](x, y), 

Ai(x,y)^Ai(x,y) + —Vir]{x,y), (7) 

where rj(x,y) is an arbitrary gauge function, defined in the whole region R. As a result, 
the variational derivatives with respect to <fi n , and A x , A y are related through fundamental 
functional identities 

= — V V; . 

5<f) n (x,y) 2e^ y SAi(x,y) 

The occurrence of such identities is typical of gauge theories.^ Moreover, identities relating 
variational derivatives appear already in some problems of classical variational calculus with 
degenerate (i. e., invariant under symmetry transformations) functionals.0 As in degenerate 
theories the number of variables exceeds the number of independent Euler-Lagrange equa- 
tions, complementary relations should be normally imposed to eliminate irrelevant degrees 
of freedom and close the system mathematically. Whereas in bulk superconductors and 
single junctions the elimination of unphysical degrees of freedom amounts merely to an ap- 
propriate choice of gauge, in periodic weakly coupled structures this problem has additional 
implications. Namely, in the presence of the Josephson interlayer coupling phase differ- 
ences $ n .n-i and $ n +i,n at two successive barriers are in themselves not independent, which 
means, mathematically, that we are dealing with a variational problem with constraints. 
Unfortunately, this fundamental feature was completely overlooked in previous literature. 

The variations with respect to f n are independent and can be taken first. Varying under 
the assumption of arbitrary 5f n at the boundaries, we obtain 
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(8) 



1 + C 2 (T) ]T [V 2 - [VA(:r, 2/) - 2eA(z, y)f 



i=x,y 



f n {x,y) - f*{x,y) = 0, 



(9) 



(x,y) e R Sn ; 



dfn 
dy 



X, Lyl) 



dfn 

dy 



[x,L y2 ) = 0; 



(10) 



— ^ (-a/2 + np, y) = — [f n (-a/2 + np, y) - / n _i (a/2 + (n- l)p, y) cos $n,n-i(2/)] , (U) 
ax z^o 



df n ot 

— - (a/2 + np, y) = -— [/„ (a/2 + np, y) - f n+1 (-a/2 + (n + l)p, y) cos $ n+ i,„(y)] . (12) 
ox z^o 

Here, equation (|) is the usual GL equation for the bulk order parameter. Relations ( [TQ|) 
are the usual GL boundary conditions at the superconductor /vacuum interfaces. Boundary 
conditions (|H]), (0), describing the suppression of the order parameter due to the Joseph- 
son currents at the superconductor /insulator interfaces, are of the type first derived by de 
Gennes0 for a single junction. 

By contrast to f n , the variables A x , A y are defined in the whole region R. For these vari- 
ables, continuity up to the second-order partial derivatives at the superconductor/insulator 
interfaces should be assumed. The corresponding Euler-Lagrange equations are 



dh{x, y) = 1 ft(x,y) 
dx 2e A 2 (T) 



d(j) n (x,y) 
dy 



2eA y (x,y) 



4irj ny {x,y), 



(13) 



dh(x,y) = 1 fn(x,y) 
dy 2e A 2 (T) 



d<f> n (x,y) 
dx 



2eA x (x,y) 



^jnx(x,y), (x,y) e R Sn ; (14) 



dh(x,y) 
dy 



4vrj /n {-a/2 + np,y) f n -i (a/2 + (n - l)p, y) sin$ n>n _i(y) = 47rj n , n _i(y), (15) 



Jo 



7C(3)ax(W0 
6 



eiV(0)eoA 2 (T), 



(16) 
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dh(x, y) 
dx 



0, (x, y) e R K 



(17) 



Here, equations (jjjD, QT^) are the Maxwell equations in the S-layers, with j n being the 



intralayer supercurrent densities. Equations fll5|), (17) are the Maxwell equations in the 
barrier regions, with j n>n -.i being the Josephson current density between the nth and the 
(n — l)th layers. Relation (|I1J) is the definition of the Josephson critical current density in 
a single SIS junction.0 

Equations (|T5|)-(|r?j) should be complemented by boundary conditions at the outer inter- 
faces y = L yl L y2 . [When deriving these equations, we have only assumed 5A x (x, L y i) = 
5A x (x, L y2 ) = 0.] As we do not consider here externally applied currents in the y direction, 
the first set of boundary conditions follows from the requirement [j ny \ y=L 1 L 2 = 0: 

d(p n (x,y) 



dy 



2eAy(x,y) 



0. 



(18) 



- y=L y i,L V 2 

Applied to Eq. (|H|), these boundary conditions show that the local magnetic field at the 
outer interfaces is independent of the coordinate x: h(x,L y i) = h(L yl ), h(x,L y2 ) = h(L y2 ). 
The boundary conditions imposed on h should be compatible with Ampere's law h(L y2 ) — 



h(L y i) = Airl obtained by integration of Eqs. (0), (0) over y, where 



Ly2 Ly2 

1 = J d yJnx(x,y) = J dyj n ^ x {y) 

Lyl Lyl 



(19) 



is the total current in the x direction. Throughout this paper, depending on a physical 
situation under consideration, we will employ three types of boundary conditions on h: 

h(L yl ) = h(L y2 ) = H, (i); 
h(L yl ) = H - 2ttJ, h(L y2 ) =H + 2ttI, (ii); (20) 
h(L yl ) = H, h(L y2 ) = H + 4tt/, (iii). 
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As usual, the Maxwell equations ([13|), (|I~4D yield the current- continuity equations inside 
the S-layers: 



Z) V * fn( x ' V) Ni<f>n(x, y) - 2eAi(x, y)\ 



0. 



(21) 



The conservation of Josephson interlayer current is readily verified from ([15]). Using Eqs. 

(|H|), ([15|) and assumed continuity of dh/dy, we arrive at the boundary conditions 

d(f) n (x,y) 



dx 



2eA x (x,y) 



fn(x,y) 



a 



- x=—a/2+np 



—fn-i (a/2 + (n- l)p,y) sin $ n , n _i (y), 

(22) 



^n(x,y) 
dx 



2eA x (x,y) 



fn(x,y) 



- x=a/2+np 



a 

Wo 



f n+1 (-a/2 + (n + l)p,y)sm$ n+lin (y), 



(23) 



reflecting the continuity of the x component of the supercurrent at the internal interfaces 
x = ±a/2 + np. 

Integrating Eqs. (0), ( pT|) over x and applying boundary conditions ([11]), flT2|) and (p2|), 
(E3|), respectively, we obtain very useful integro-differential representations 



U(x, y) - mx, y) - C 2 (T) £ [V,0 n (x, y) - 2e^(x, y)] 2 / n (x, y) + C 2 (T) 



2/rTA d 2 f n (x,y) 



i=x,y 



dy 2 



«C 2 (T) 
2< 



[/„ (-a/2 + np, y) + /„ (a/2 + np, y) 



-fn+i (-a/2 + (n + y) cos $ n+ i >n (y) - / n _i (a/2 + (n - y) cos $ njTl _i(y)] , (24) 
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£(*,v) |^f^-2e^(x,y) 



a 



2a£ 



<9y 



[/n (- fl / 2 + "P, y) /n-i (a/ 2 + (n - y) sin $ n , n _i(y) 



-/n ( a / 2 + n P> y) fn+i (-a/2 + (n + l)p, y) sin $ n+1>n (y)] 



(25) 
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a/2+np 

where (x, . . .) = - / (x, . . .) dx denotes averaging over the interval —a/2 + np < x < 

—a/2+np 

a/2 + np. 

By summing Eqs. (]25 ) over the layer index n, integrating and applying boundary con- 



ditions (|18D , we obtain the integral 



ti( x >y) 



d<t> n (x,y) 
dy 



2eAy(x,y) 



0. 



(26) 



which is, physically, the conservation law for the total supercurrent in the y direction. 
Mathematically, equation (|26| ) has the form of a constraint relation between variables dcj) n / dy 



and y4 y £3c3 To find the rest of constraints of the theory, closing the system of equations, we 
must minimize the functional (|5|) with respect to <f> n and d<p n /dy. 

By virtue of fundamental identities (|f), a naive variation of @ with respect to (f) n (with 
arbitrary 5(fi n at the boundaries) does not yield new equations. Indeed, the correspond- 
ing Euler-Lagrange equation reduces to the conservation law (^Tj), while surface variations 
merely reproduce boundary conditions (|T8|), and (f22|), (f23|). Considering variations of the 
type <f) n {x,y) -> (f) n (x, y) + eip n (y), d<f) n (x,y)/dy -> d<f) n (x, y)/dy + edijj n {y)/dy, where e is a 
small parameter and ip n {y) are arbitrary functions of y, we arrive at Eqs. (|25|). To obtain 
genuinely new equations, minimizing (||) with respect to (fi n and dcf) n /dy, we must enlarge 
the class of allowed variations. 

A mathematically rigorous approach to this problem is as follows. While varying 
<fr n (x,y) -> (j) n (x,y) + 5(f) n (y), d(p n (x,y) /dy -> d(j) n (x, y)/dy + d5<f) n (y)/dy, with 5<f) n (y) 
being small arbitrary functions of y, instead of integrating by parts, we impose additional 
constraints 



d<t> n (x,y) 

dy 



2eA y (x,y) 



(27) 
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compatible with boundary conditions (|18|) and constraint relation (p6|). The requirement 
of compatibility with the current-conservation law fl2"5|) automatically yields another set of 
constraints 

f n (-a/2 + np,y) f n -i (a/2 + (n - l)p,y) sin $ ni „_i(?/) 

= /„ (a/2 + rip, ?/) f ri+l (-a/2 + (n + l)p, y) sin $ n+ i,„(y). (28) 

The above procedure is formally equivalent to minimization of @ with respect to indepen- 
dent variations of <fi n and dcj) n /dy. As this class of variations of <p n and d<p n /dy is larger than 
that employed in deriving fl25|), we can argue that Eqs. (^7|) and (^) provide the sought 
necessary conditions for the true minimum of the free-energy functional ([5]). 

The physical meaning of Eqs. ([27]) and (|28| ) is quite transparent. Constraints (^) mini- 
mize the kinetic-energy term in ([5]) with respect to variations d(j) n (x, y)/dy — > d(p n (x, y)/dy+ 
5ip n (y), where 5ip n (y) are small arbitrary functions of y. They show that the average in- 
tralayer currents in the y direction are always equal to zero, and, as a result, 

h(-a/2 + np,y) = h(a/2 + np 1 y) . (29) 

[See Eq. fljjP ]- These constraints appear already in the case of decoupled S-layers. 
By contrast, constraints (^) are uniquely imposed by the Josephson interlayer coupling. 
Their function is to make the Josephson energy stationary with respect to variations 
4> n (x,y) — > 4> n (x,y) + S(p n (y) and to assure the conservation of the total Josephson cur- 
rent / in neighboring barriers [see Eqs. (0), (P^QI- 

As no other conditions are imposed on the variables, we can satisfy ( |2"8"D by choosing 

f n (x,y) = f n -i(x-p,y) = f(x,y), f(x + np,y) = f(x,y), (30) 
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$n+l,n(y) = $n,n-l(y) = Hv)- (31) 

These relations finalize the determination of a closed, complete, self-consistent system of 
mean-field equations for a S/I superlattice in the GL regime. 

Constraints Q27D, ( PBQ and their corollaries (^)-(BT|) belong to key results of this paper. 
Derived by means of a rigorous mathematical analysis of the impact of gauge invariance, 
they are not restricted to the functional @, but should hold for any superconducting weakly 
coupled periodic structure. To illustrate their importance, we point out that Eqs. (|28|) , (p9|). 
for example, completely rule out any possibility of single Josephson vortex penetrationi'i and 
triangular Josephson vortex lattice,! proposed without appropriate physical and mathemati- 
cal justification. On the contrary, they imply that the distribution of the local magnetic field 
due to the Josephson vortices has, in general, the periodicity of the multilayer, as recently 
verified experimentally.^ It should be noted, however, that although the role of constraints 
(fZ5|), fl2P| ) in minimizing the free energy and closing the system of Euler-Lagrange equations 
for /„, <j) n , and A has not been realized until now, relations fl3U|), (PT|) were implicitly em- 
ployed in phenomeno logical calculations of i? C 2oo- Moreover, relations of the type 
(|7|), (|30|) , ( |3TD were used by TheodorakisS in his particular exact solution of the LD model 
in a parallel field. 

The equations of this subsection admit exact solutions in two limiting situations: the 
single-junction case, when a ^> max{£(T), A(T)}; the thin-layer limit, when £ <C a <C 
min{C(T), \(T),a~%, W} (W = W y = L y2 - L yl is the length of the S-layer in the y 
direction). The single-junction case is well-known. The thin-layer limit will be extensively 
discussed in the next subsection and in Sec. III. 
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B. The thin-layer limit 

The mean-field equations of the previous subsection allow remarkable simplification in 
the thin-layer limit, when £ *C a <C min{£(T), A(T), a _1 £o, W}. 

First, we can neglect the x-dependence of /, defined by Eq. fl5D|): f(x,y) = f(y). Second, 
fixing the gauge by the condition 

A x (x, y) = 0, Ay(x, y) = A(x, y), (32) 



we can neglect the x-dependence of 0„ as well: <p n (x,y) = 4> n {y)- In the gauge fl32|), 
^n.n— 1 (y) = 0n(y) - <Pn-i{y), and Eqs. (|3Tj) become 

<fin+i(y) + <frn-i(y) = 2<frn(y), <frn(y) - <frn~i(y) = <f>{y), 

with the solution 

<f>n(v) = ntf>(y) + x(y), (33) 

where 4>(y) is the coherent phase difference (the same for all the barriers), and x{v) is an 
arbitrary gauge function, allowed by particular gauge transformations 4> n (y) — > (f> n (v)+x(y)j 
A(x,y) - A{x,y) + i&M Without any loss of generality, we can set x = 0. 

In view of independence of / and n from x in the thin-layer limit, the physical meaning 



of constraints (p7|), (|28|) becomes even more obvious. Thus, Eqs. (p8|) are now the conditions 
of stationarity of the Josephson energy with respect to all allowed variations of <p n . Due to 
Eqs. fl2~T|), the term in ( p4| ) responsible for the kinetic energy of the intralayer currents 
becomes 



C 2 (T) £ [VA(x, y) - 2e^(x, y)f /(x, y) - ( 2 (T) 



i=x,y 



-2eA(x,y)l /(y) 



cfy 
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C 2 CO 



d(j)n{y) 

dy 



-2eA{x,y) 



f(y)+4e 2 ( 2 (T) A*(x,y) - A(x,y) f(y) 



4e 2 ( 2 (T) A 2 (x,y)-A(x,y) f(y) 



which shows that conditions (P7|) minimize the kinetic energy for a given configuration of 
the vector potential A. 

Concerning the Maxwell equations, the right-hand side of Eq. ([13|) is of order a 2 /A 2 (T) 



and can be discarded. Equation ([14]) can be altogether dropped. Thus, we arrive at a closed 
set of equations 



d 2 A(x, y) , . „ 
4^ = 0, (x,y)eR; 



dx 2 



1 d 2 A(x,y) 
An dydx 



jof (y) sin 4>{y), (x,y) e R bn ; 



dy 



dA(x,y) 
h{x,y) = — — — , [x,y) G R., 



(34) 



(35) 



(36) 



(37) 



f(y)~ f(y)-4e 2 ( 2 (T) A*(x,y) - A(x,y) f(y) + C(T) 



d 2 f(y) 
dy 2 



«C 2 (T) 



[1 - cos <f){y)]f(y), (x,y)ER Sn . 



(38) 



These equations, of cause, should be complemented by continuity conditions on A, dA/dx 
and boundary conditions ( |T0l) and (|20|). 

It is worth noting that an immediate consequence of Eqs. (|34|), (|38]) is independence 
of the local field /i from the coordinate x in the whole region R: h(x,y) = h(y), — oo < 
a; < +oo. This result is fully compatible with the requirement (^) and demonstrates that 
the intralayer supercurrents in the thin-layer limit are unable to screen out the magnetic 
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field: The situation is very familiar from the physics of isolated superconducting films with 
a<A(T)£0§ 

Our next objective is to eliminate the vector potential and obtain a closed set of equations 
involving only / and <fi. Equations Q25D, fl3"5] ) can be easily solved for A in the nth "elementary 
cell" R n = R Sn U Rb n (the S-layer plus the adjacent barrier). Applying the continuity 
conditions on A, dA/dx, boundary conditions fl2Cp and the constraint relation Q36D , we get 



A(x,y) 



47TJ0 J duf 2 (u) siacj)(u) + Hi 

Lyl 



n d(f>{y) 

( x -n P )+-—, (39) 



where H X = H for (gDJ), (i), (iii), and H X = H- 2nl for (^DJ) , (ii). Matching Eq. Q^) to an 
analogous solution in the adjacent cell R n -\ leads to the solvability condition 

= 8nej p [ duf 2 (u) sin0(w) + 2epH 1 . (40) 
rfy J 

Equation ( ^0|) is nothing but an analogue of the Ferrell-Prangei relation for a single Joseph- 
son junction, which can be readily verified by differentiation. From this point of view, the 
quantity (87rejop) -1 ^ 2 should be identified with the Josephson penetration depth Aj. Note 
that instead of the factor 2A entering the definition of Xj in the single-junction case,0'ii in 
our case we get the period p. 

With the help of Eq. (|3^), we arrive at the expression for the vector potential in the 
whole region R = U+^^i^ = R s U R b : 

M*,V) = ^^-*, (x,y)eR. (41) 



This equation should be substituted into Eqs. (|37|), (|38|). 

In this manner, we obtain a closed, complete set of equations describing a thin-layer S/I 
superlattice in an external parallel magnetic field: 



A(x, y) = A 00 f(y)J2S Rsn (x, y) exp [intf>(y)] 



(42) 



$R Sn (x, y) 



1, for (x,y) e R Sn , 



0, for (x,y) g R Sn ; 



P 

c< 2 (T) 
< 



d<P{y) 

dy 



f(y) + C 2 m 



d 2 f{y) 
dy 2 



f{y) 



[l-cos^(y)]/(j/) = 0, 



(43) 



(44) 



(45) 



Aj = (8irej p) 1/2 



(46) 



i rf0(y) 

2ep oh/ 



(47) 



= j x (y) = jo f 2 (y) sin 0(2/) 



1 d%) 

47T (iy 



!/-' 



(48) 



with boundary conditions ( p0|) and I = / dyj(y), where j(y) is the x component of the 



supercurrent density (both in the S-layers and the barriers). The y component of the in- 
tralayer supercurrent, whose average over the layer thickness is equal to zero, within the 
accepted accuracy enters the theory only implicitly, via the average kinetic-energy term in 
Eq. ©. 

Significantly, the coherent phase difference (the same for all the barriers) obeys only one 
nonlinear second-order differential equation flj5p with only one length scale, the Josephson 
penetration depth Aj, as in the case of a single junction.ll'S Due to the factor f 2 , equation 



19 



(f45p is coupled to nonlinear second-order differential equation (f43|), describing the spatial 
dependence of the superconducting order parameter / (the same for all the S-layers). In 
the latter equation, the term proportional to a 2 /p 2 accounts for the average kinetic energy 
of the intralayer currents, while the term proportional to a accounts for the kinetic energy 
of the interlayer Josephson currents. The Maxwell equations P7|), (pT^D, combined together, 
yield Eq. (|45|), as they should by virtue of self-consistency. 

It is instructive to compare the above equations with those now circulating in literature 
concerned with the phenomenological LD model. As already mentioned, neither mutual 
dependence of the Euler-Lagrange equations for <ft n and A, nor fundamental complemen- 
tary relations of the type fl2T|), (28), minimizing the free energy, have been established 



in previous publications. Left with an incomplete set of equations, some authors make 
a non-self-consistent approximation /„ = 1 and, regarding the phase differences, propose 
a mathematically ill-defined infinite set of differential equations with two different length 
scales. [See, e. g., Refs. (ifl).] In view of the conditions (PT|), these equations reduce to our 
Eq. (H) with f = 1. 

Finally, the free-energy functional @ in the thin-layer limit after a transition to the 
mean-field approximation with respect to A takes the form 



n[f,(f>;H} 



An 



j »2 



W X W Z / dy 



dm 

dy 



C 2 (T) 


0' 


d(f)(y) 
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dy 



f\y) + 



aC\T) 



+4e 2 C(T)X\T) 



1 d<j>(y) 
2ep dy 
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[l-coscKi/)]/ 2 ^) 



H 



(49) 



where W x = L x2 — L x \. As expected, minimizing Eq. (|49|) with respect to / and the phase 
difference 0, and neglecting terms of order a 2 /X 2 , we arrive at Eqs. (f43|)-(f45|). 



The functional ([EJ) and complementing Maxwell equations (f48|), (pi|) contain much more 
physical information than the phenomenological LD model in a parallel field: the domain 
of validity is exactly determined, all the coefficients are microscopically defined, and a finite 
S-layer thickness is explicitly taken into account. (As we show in Sec. Ill, this factor removes 
unphysical divergence of i? C 2oo, typical^! of the LD model.) Another important difference is 
the proportionality of the condensation energy in fl49|) to the layer thickness a, instead of 
the period p in the LD functional. 

The equations of the thin-layer limit admit exact solutions for all physical situations of 
interest. These solutions are the subject of the next section. 

III. MAJOR PHYSICAL EFFECTS IN THE THIN-LAYER LIMIT 

A. The Meissner state in a semi-infinite multilayer. The superheating field 

H s = (epXj)- 1 

Consider a semi-infinite (in the y direction) multilayer with L y \ = 0, L y2 = +00 in the 
external fields 

0<H<H s = (ep\j)-\ (50) 

with boundary conditions of the type (pD|), (iii): 

h(0) = -L^(0) = #, Z^+oo) = -L^(+oo) = H + 4tt/ = 0, 0(+oo) = 0. (51) 
2ep ay 2ep ay 

For 
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«C 2 (T) 



«1, 



(52) 



arc 



the Meissner solutions of Eqs. ([4q)-(f 

4>(y) = — 4arctan 



E exp 



JL 

Aj 



h s + je 2 - e 2 



%) 



2M a 



# s 2 - i7 2 



cxp 



JL 



H s + ^E 2 - E 2 \ + E 2 exp 
EE. 



2y 
A./ 



(53) 



(54) 



j(y) 



2vrA„ 



E 2 S - E 2 



[E s + yj E 2 - E 2 


— E A exp 


2y~ 

. Aj. 


exp 


y 

. Aj. 




E s + JH* - E 2 


2 r 

+ E 2 exp 


2?/" 
A./. 


2 



(55) 



f(y) = i - 



4«C 2 (T) H 2 


E s + y/H* - E 2 


r 

exp 


22/" 
Aj. 




a£ 


E s + ^E 2 - E 2 


2 

+ E 2 exp 


. Aj. 


1 2 



(56) 



The Meissner solutions persist up to the field H s = (epAj) -1 that should be regarded as the 
superheating field of the Meissner state. 

Indeed, as we will show below, the presence of Josephson vortices inside an infinite 
multilayer becomes energetically favorable at a field H = Ed^ < E s . As in the case of 
the well-known Bean-Livingston bamerfflS in semi-infinite type-II superconductors, the 
penetration of Josephson vortices at fields i? c ioo < E < E s is prevented by a surface barrier 
due to surface currents j(0). [Compare the discussion of the superheating field in the case 
of a singe junction,il where it is given by the expression H s = (2eAAj) -1 ]. Equation fl55|) 
shows that |j(0)| increases in the interval < H < H s /y/2, reaches its maximum value 
at E = H s /y/2, decreases in the interval E s j\f2 < E < E s and vanishes at H — E s . 
Moreover, the phase difference at the surface 0(0), being a non-positive, monotonously 
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decreasing function of H in the whole interval < H < H S} also reaches its minimum value 
0(0) = — 7r at H = H s . The appearance of the phase difference -it can be attributed to the 
formation of a line singularity of the amplitude of condensation (ifj^(r)ipi(r)) at the outer 
interface of the barrier ("the Josephson vortex core"). In addition, the magnetic flux per 
"elementary cell" at H = H s is $ =$o/2, where $o = 71 l e is the flux quantum. 

Finally, from the second of Eqs. (]5lf) and the condition H < H s for the Meissner 
solutions, the maximal value of the total Josephson current |J| = J max in a semi- infinite 
multilayer is 

/max = ~7~~ — 2\jj . (57) 
47T 

For |/| > J max , the field at the boundary is h(0) > H s , and the stationary flow of the 
Josephson current is disrupted by the penetration of Josephson vortices that move under 
the influence of the Lorentz force. 

Thus, in fields H > H s , only vortex solutions are possible. Owing to the specific feature 
of the thin-layer limit, i. e., the absence of screening by the intralayer currents, the "tails" of 
magnetic field distribution of individual Josephson vortices overlap in the layering direction, 
causing the formation of unique vortex structures that we term here "Josephson vortex 
planes". We begin the discussion of these structures form a single "vortex plane", forming 
in an infinite layered superconductor at the lower critical field H = i? c ioo. 

B. The lower critical field i^doo = 2(7reAjp) _1 in infinite multilayers. Vortex planes 

Consider now an infinite (in the y direction) layered superconductor with L y \ = — oo, 
L y 2 = +oo, subject to boundary conditions of the type (|2"0"D, (i), with H = 0. The condition 
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(|52|) is supposed to be fulfilled. We are looking for a vortex solution with one flux quantum 
$o P er "elementary cell", i. e., with 



f+oo) — d>(— oo) = 2n 



1 



2ep dy 



(±00) = 0, 0(0) = n. 



(58) 



The sought solution has the form of a kink: 



4>{y) = 4arctanexp 



y_ 



(59) 



This solution describes a single vortex plane positioned at y — 0. [Compare the phase 



difference 0(0) = n of Eq. (|59p with the phase difference 0(0) = — 7r of Eq. (^) at the 
surface of a semi- infinite superconductor in the field H = H s , when a vortex plane only 
starts to penetrate. After the actual penetration, the phase difference changes by 2tt, as 
expected from general considerations.il] 

Corresponding distribution of the local magnetic field is given by 



h(y) = (epXj) 1 cosh 1 



Notice that at the vortex plane h(0) = H s . The density of the Josephson currents is 



(60) 



j{y) = -2j cosh 



-2 



y_ 

Aj 



sinh 



y_ 

Aj 



(61) 



At the vortex plane, j(0) = 0. The Josephson currents vanish exponentially at y —>■ ±00 and 
reach their peak values at y = ± ln(l + V2)Xj ~ ±0.88Aj. As regards the order parameter, 



we get 



f(y) = 1 - 



4a( 2 {T) exp 



2y 

Aj 



a£ 



1 + exp 



2y 



2 ' 



(62) 



Notice that Eqs. (|5"DD-(|5"2"D, considered in the half-space < y < +00, have exactly 



the same form as the solutions (p4|)-(pq) for a semi-infinite multilayer in the external field 
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H = H s , in full agreement with our interpretation of H s as the penetration field for a single 
vortex plane. 

To find the lower critical field at which the solution (|59D becomes energetically 

favorable, we must consider the free-energy functional (P?D, which in this case takes the 
form 



n[<j>(y);lT\-n[H] Nv=0 = N oM W z 



Jo 
2e 



dy 



1 — COS0(?/) 



A 2 , 



dy 



1 [(f>{+oo)-<f)(-oa)}H 



4tt 



2e 



where [//] 



(63) 



AT„=0 



0; H] is the free energy in the absence of vortices (N v is the 



number of vortex planes), and N ce u = W x /p is the number of "elementary cells". Inserting 



into (63), we obtain the free-energy contribution of a single vortex plane: 



Q\H) 



N v =l 



tt[H} 



Nv=0 = N ceU W z 



4Ajj _ $ H 
e 4tt 



(64) 



with $ — 7r/ e the flux quantum. From the condition Q [if cloo ] 7V ^ =1 
lower critical field is 



coJJV«=0> 



the 



-f^cioo = 2(7iep\j) 



-i 



2 $ r 



7T 7TJ)A„ 



(65) 



as in the case of a single junction, apart from the factor p in the denominator instead 
of 2A(T)ffl As expected, H cloo = 2H s /tx < H s = h(0). On the contrary, Eq. (H) is 
completely different from previously proposed ones for layered superconductors,! based on 
an invalid assumption of single- vortex penetration. 

From the proportionality of the right-hand side of fl6"4|) to N ce u, we infer that the total 
number of Josephson vortices (i. e., ID singularities of the amplitude of condensation) in one 
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vortex plane is equal to the total number of elementary cells. This means that Josephson 
vortices penetrate all the cells simultaneously and coherently. As in the case of a single 
junction,il the quantity 

E = ^ (66) 

can be identified with the self-energy of a single Josephson vortex per unit length (in the z 
direction) . 

In higher external fields (H ^> // c ioo), we expect to get a "stack" of N v vortex planes 
with the total number of Josephson vortices N Vtot = N v N ce u. (See Fig. 2.) 

C. The vortex state in intermediate fields -ff c ioo <C H <C [ea£(T)]~ . The lower critical 
field H c i\y = Ti/epW in finite-size samples (W <C Aj). A paramagnetic effect 

Consider a finite-size (in the y direction) multilayer with — L y \ = L y 2 = L, W = 2L, in 
the field range i? c ioo *C H [eaC(T)] -1 and in the absence of externally applied current 
(1 = 0), i. e., subject to the boundary conditions (|20|) , (i). The validity of the condition 
(152]) is again assumed. [The upper bound [ea£(T)] _1 for the field range means that we are 
concerned with H -C H c2oo (T).] 

Under these assumptions, the phase difference up to first order in the small parameter 
(ep\jH)~ 2 is 

f-iF 

(j>{y) = 2epHy + ttN v (H) - -f / [sin (2epHy) - 2epHy cos (epWH)} . (67) 

A(epAjH) z 

The constant of integration ttN v (H) accounts here for the phase shift due to N v vortex 
planes [ir per each vortex plane, see the last of Eqs. (§§)]. The number of vortex planes N v 
is itself a singular function of the applied field H: 
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N V (H) 



\epWH- 






7T 







(68) 



where [u] means the integer part of u, and $ = pWH is the flux through an "elementary 
cell" . This choice of the constant of integration guarantees that the energy of the Josephson 
coupling Ej in takes the minimal value for a given H: 



EAH]= ?m w , ww °?v) 



1 $ 


7T$ 






sin — — 




$0 


$ 





Aix a£ 

(This expression should be compared with its analog for a single Josephson junction, i) 
The physical quantities corresponding to ( |67|) are 



(69) 



cos (2epHy) — cos (epWJT)] 



f(y) = i - 



«C 2 (T) 
2< 



A(epXjH) 2 
j{y) = (-1)^0 sin (2ep#y) 
-l) % cos(2epify) 



(70) 



(71) 



V2epHQ{T) |sin (epFW) | cosh ^ 



l + 2[eptfC(T)r 



1 + 2 [eptfCCO] sinh 



w 



(72) 



[The term 2 [epifC(T)] 2 in the denominators of ( 715 ) can only be retained if p ^> a. 
In the limit W > ((T), \y\ < W/ 2 > equation Q72|) becomes 



/(y) = i 



«C 2 (T) 
2af 



cos (2ep/fy) 



(73) 



l + 2[epHC(T)Y 

Equations of the type © and fl73|) for N v = 2m (m is an integer) were first obtained 

by TheodorakisS in the framework of the LD model. Our equation fl67D for N v = 2m should 
also be compare with an analogous solution for an infinite single junction given, for instance, 
in Ref. @). 
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The singular function N V (H) introduces discontinuities in Eqs. (|67D, ([70|)-( [f^) . These 
discontinuities witness that the system undergoes a first-order phase transition when a vortex 
plane penetrates or leaves the sample. [Compare with the discussion of a single junction in 
Ref. ^ 



The positions of vortex planes y v correspond to local maxima of the field h(y) in Eq. 
(fTOf). [In the case (f73|), y v exactly coincide with local minima of f{y)-\ In the vortex planes 
y = y v , the microscopic magnetic field is higher than the applied one: 



h(y v ) = H 



1 + 



4(epXjH) 2 



1 + (-l)^cos (epWH) 



(74) 



which is expected for any vortex solution. The Josephson current density j{y) = 
vanishes both in the vortex planes y = y v and in the planes of local minima of h(y), y = 
y v ±n/2epH. When passing through zero in these planes, j(y) changes the sign, as depicted 
in Fig. 2. 

From Eq. ( |B"BD with N V (H) = 1, we obtain the lower critical field H c i\y in a finite 
multilayer with W <C Aj: 

vr it 2 „ \„ 



H. 



clW 



epW~ 2^ cloo ^ >>Fclo °- 



(75) 



The definition of the magnetization M,0 

1 +L 

AttM = — J dyh{y) - H 



(76) 



and Eq. (70) yield 
M(H) -- 



1 



16nH(ep\ 



I sin (epWH)\ 



-l) Nv cos (epWH) 



(77) 



epWH 

The magnetization (|77J) shows distinctive oscillatory behavior and discontinuities at 
epWH — ► ttN (N is an integer), when a vortex plane penetrates or leaves the sample. 
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Interestingly enough, the right-hand side of (ff7|) passes through zero and may have both 
signs. Thus, for $ = pWH ^> $ , the sample exhibits a small paramagnetic effect, if 



N v $ <$ < (N v + \ - $ 



M(H) 



l6nH(ep\j) 2 



| cos (epWH)\ 



| sin (epWH)\ 
epWH 



> 0. 



(78) 



D. Fraunhofer oscillations of the Josephson current in multilayers with W -C Aj, in 
the field range < H <C [ea£(T)]~ . "Edge pinning" of the vortex planes 



Now we proceed to the case of a finite-size (along the layers) multilayer with —L 



2/1 - 

L y2 = L, W = 2L in the presence of an externally applied current /, i. e., subject to the 
boundary conditions (|20|) , (ii). (Compare the discussion by Owen and Scalapincll of the 
single-junction case.) The relation fl52[) is supposed to hold. The applied magnetic fields are 
within the range < H <C [eaC(T)]~ . 

Assuming W <C Aj, we can consider W 2 /X 2 as a small expansion parameter in Eq. (|45|). 
In this way, we obtain 

0(y) = 2epHy + 7rN v (H)+(p 





W^ 2 




2 p 


4 


Aj 







7T$ 

sin (2epHy + <p) — 2epHy cos — — cos 99 — sin </? 

^0 



(79) 



sin 



7T<3> 
$ 



(80) 



%) = # 



W 2 



4 Aj 



"$0~ 


2 r 







7T$ 

cos (2epHy + cp) — cos — — cos <p 



f(y) = 1 - 



«C 2 (T) 
a£ 



1 - 



(-l) Nv cos(2epHy + (p) 
l + 2[epFC(T)] 2 



(81) 
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V2epH((T) |sin (epHW)\ cos <p cosh g 

~ m n T~\ 



1 + 2 [ep#CCO] sinh 



x/2C(T) 



jj^yg^ffCg) cos (epHW) sin y sinh ^ 



(82) 



l + 2[epHC(T)Y cosh^^y 

where N V (H) = SEWJL [ s the number of vortex planes, $ = pWH is the flux through an 
"elementary cell", $o = ^/ e i as usual, and the constant tp (\tp\ < 7r/2) parameterizes the 
total Josephson current / given by (|80|). Equation fl8"0|) yields the well-known Fraunhofer 
pattern, the only difference from the single-junction case being the occurrence of the period 
p in place of 2A(T).il@ Note that the first zero of the Fraunhofer pattern, by virtue of ([75]), 
corresponds to the lower critical field H c iw- In the absence of the transport current, i. e., 
for ip = 0, equations (|7^) , (|Sll), ( [32"1) reduce, respectively, to (|67|), (f70|) and (|72|) , as they 
should. 

The self-consistency of our calculations can be easily verified by means of Ampere's law 
h(+L) — h(—L) = 4irl. It is assured by terms proportional to W 2 /\ 2 j in (|79|) and (jSTp that 



explicitly take into account the effect of self-induced fields. Although Eq. (|80D was first 
derived in the framework of the LD model in Ref. (i), the authors of that publication, based 
on an incomplete set of equations, were unable to find the phase differences self-consistently 
and evaluate the local magnetic field in first order in W 2 /Xj. As a result, they arrived at 
an absolutely incorrect conclusion that Fraunhofer oscillations of I could be observed in the 
absence of Josephson vortices. Unfortunately, this misunderstanding is shared in some other 
recent publications.^ Therefore, we provide below a detailed and rigorous clarificatk 



cion. 



As we see from Eq. (|8lD , in the presence of the transport current /, the vortex planes 



are shifted by the Lorentz force to new equilibrium positions [local maxima of h(y)] : 
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Vv = Vv 



2epH' 



(83) 



where y v correspond to local maxima of the right-hand side of (]81|) for ip = 0. The local 
magnetic field in the vortex planes now is 



h(y v ) = H 



1 + 



1 



A{ep\jHy 



1 + 



\N V 



cos (epWH) cos ip 



> H. 



(84) 



In equilibrium, the Lorentz force fi per elementary cell acting on the vortex planes is 
counterbalanced by the pinning force f pin that can be defined as0 



fpin\Y) 



1 dU ptn (Y) 
N cell W z dY 



(85) 



where U P i n (Y) is the pinning potential arising owing to the shift by Y of the vortex planes 
from their equilibrium positions in the absence of the transport current I. To evaluate the 
pinning potential, we consider the increase of the free energy in first order in a( 2 (T)/a!; , 
caused by such a shift. Noting that first-order corrections to f(y) ~ 1 and h(y) w H do not 
contribute to the free energy, taking <p(y) w 2epHy + nN v (H), making the transformation 
y — > y — Y and substituting into Eq. (f03), we obtain 



U pm (Y)=n[H;Y]-n[H;Y = 0] 



N C£ll WW; 



2e7r$ 



sin 



7T$ 



— COS f 



2tt$ F 



(86) 



It is very instructive to rewrite Eq. (pq) as 



1 $o 



f/ pm (F) = JV^W,-^ [2j(+L; Y = 0)+ j(-L; Y) - j(+L; Y)] 



(87) 



where 



j(±^; Y) = (-l) Nv j sm(±2epHL + 2epHY) 
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are the surface currents in the presence of the shift Y. We see that the pinning potential for 
\Y\ < ir/4epH arises owing to the emergence of additional surface currents on the opposite 
sides of the superconductor. At 2epHL = nN + 1 (N is an integer), j(—L; Y) = — Y), 
i. e. these currents flow in the opposite directions, and the pinning potential reaches its 
maximum. On the contrary, at 2epHL = ttN (N is an integer), j(—L;Y) = j(+L;Y), 
i. e. the surface currents flow in the same direction and mutually compensate each other 
in (|87D , the pinning potential vanishes, and vortex planes freely penetrate or leave the 
sample. (Compare with the discussion at the beginning of this Section of the case of a 
semi-infinite multilayer for H = H s .) The surface currents also flow in the same direction 
and mutually compensate each other when the magnitude of the shift \Y\ reaches the value 
\Y\ = Ymax = 7f/4epH. Moreover, the pinning potential vanishes for $ 3> $o- 
From ([85]) and (|86|), we obtain the pinning force for the shift Y: 

/27T$ Y \ 

W^-'trT' (88) 



2tt$ Y 
o 

7T<I> 

sin — - 

$ 



sm 



2vr$ Y 



From these expressions we infer that the maximal value of the pinning force \f P i n \ for given 
flux $ is f™*=\l (|f| ;$)|$. 

In the presence of the transport current !(<£>; $) [Eq. (|80D1, the shift of the positions of 
the vortex planes, according to (p3|), is Y = —ip/2epH, with the maximal equilibrium value 
\Y\ = F miB = n/AepH. Taking into account the fact that in equilibrium = —f p i n , we 
arrive at the expression for the corresponding Lorentz force: 

& = (89) 

This expression was to be expected from general considerations,0 which prescribe for the 
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magnitude of the Lorentz force the relation|/i| = |J| where I is the transport current. It 
is therefore absolutely clear that the stationary Josephson effect becomes impossible if the 
magnitude of the transport current |/| exceeds the value / max = I ( § ; $J , because in this 
situation |/l| > \fpin\, and the vortex planes are completely depinned. 

Notice that the physics of the Fraunhofer pattern in single junctions was discussed in 
terms of a series of first-order phase transitions due to successive penetration of Josephson 
vortices long ago.@ A qualitative explanation by means of the "edge pinning" was proposed 
in the book by Tinkham.El In general, the pinning of Josephson vortices in weakly cou- 
pled superconducting structures with W <C Xj is completely analogous to the pinning of 
Abrikosov vortices by the edges of a thin [compared to A(T)] type-II superconducting film.ll 

Finally, we observe that the magnetization in the presence of the transport current 
I((p; $), according to (|76| ) and is given by 



M{H) 



16nH(ep\j) 2 



\am(epWH)\ , v 



cos (epWH) 



cos (p. (90) 



epWH 

For ip = 0, equation fl§0|) reduces to For $ = pWH > $ , N v $ <$ < 

^N v + \ — Jr^) $0; we again obtain the paramagnetic effect. [Compare with Eq. flTSP-l 

E. Critical parameters of an infinite multilayer: T coo , -£f C 2oo 

At the point of the second-order phase transition to the normal state, f 2 can be consid- 
ered as a small parameter. Thus, the term f 3 in Eq. (]68D and the right-hand side of Eq. ( |70|) 



can be dropped. Applying boundary conditions fl20D , (i) yields <p{y) = 2epHy + tcN v (H) 



With this phase difference, the linearized version of (|T^) can be transformed into 
d 2 f(t) 



dt 2 



+ 



A(T,H) - q(H) cos 2t\ f(t) = 0, (91) 
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1 - \e 2 H 2 a 2 ( 2 (T) - ™ a 

A(T, H) = 3 - —— 2 2^2-, q(H) = 

[ep((T)H} 2 2a£ (epH) 2 ' 

where we have introduced a dimensionless variable t = epHy: f(t) = f(t/epH). Hence one 
gets two independent equations: for the odd N v = 2m + 1 (m = 0, 1, 2, . . .) and the even 
N v = 2m number of vortex planes. Both of them have the usual form of Mathieu equations. 
[See Appendix A.] As to the boundary conditions, it is convenient to take — L y \ = L y2 = 
L = W/2 and, by symmetry, consider Eq. fl91| ) in the interval < y < L, with 

|(0, = |(e pM )=0. (92, 

The critical parameters T c and H c2 are now determined by the smallest eigenvalue of the 
boundary problem (0), (|92j). 

In an infinite in the y direction multilayer (L — > oo), the only bounded at the in- 
finity solutions of (|91"D are periodic Mathieu functions, with fN v =2m+i(t) occeo(t,q) and 
fN v =2m(t) occeo(7r/2 — t, q) corresponding to the smallest eigenvalues ao(q) and ao(— q) = 
a (q), respectively. Thus the critical parameters are given by the equation 

[MT,H)} coo = [a (q)] coo , (93) 

where one should fix H to obtain T coo or, alternatively, fix T to obtain H c2oo . As in the 
case of Eq. ([7^), local minima of the reduced order parameter f(t) in fliTTD correspond to 
the positions of the vortex planes: in conventional units the distance between two successive 
minima is Ay v = ir/epH, which gives the flux $ = Ay v pH = $o per single vortex. 



1. The critical temperature T c 



For magnetic fields H in the range 
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a£oX 1/2 (WO 



(94) 



the general expression for T coo resulting from Eq. (p3|) is 



T C 0C T' 



c0 



I e 2^2 2 + « / F) 2 / « 

3 < KP } °\2a^(epH) 2 , 



(95) 



In weak fields 



H < 



ep y 2a£ 7rp y 2a£ 



(96) 



pair-breaking effect of intralayer supercurrents is unimportant, and we get 



T C 00 T C Q 



7\/2C(3), 



a 



12 



<0 



(97) 



In strong fields 



vrp A/ 2a£ c 



<ox 1/2 (WO ' 



(98) 



Eq. fl9~5|) becomes 



T. 



c0 



I - 



^ 2 x (WO 



3 a^o 



(99) 



The term proportional to a takes into account pair-breaking by the Josephson currents, 
locally equal to the critical ones.0 In the absence of weak coupling (a = 0), equation fl9~9"|) 
reduces to the well-known expression for an isolated thin superconducting film.illlil 



2. The upper critical field -Hc2oo 



For a fixed T, equation fl93|) yields an implicit expression for H C 2oo as a function of T: 



[eK(mc2oo] 



3 Vp7 V2<o[eK(T)^ c2o o] 2 , 



< 2 CO 

a£ 



(100) 
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This expression exhibits the so- calledS " 3D-2D crossover" , experimentally verified, for exam- 
ple, on Nb/Ge multilayers.0 The "crossover temperature" T* can be conventionally defined 
by the relation a( 2 (T*)/a(,o = 1- 

For temperatures close to T c0 , when 



a£ 



> 1, 



(101) 



H c 2oo (T) 



12 



V2np((T) v^C(T) 7V2((^V^P^X (£o/0 



1 - 



T 



-c0 



(102) 



In this " 3D" regime, the positive kinetic energy of small interlayer Josephson currents in fl49|) 
competes with the negative intralayer condensation energy. The superconductivity of the 
S-layers is strongly depressed by the vortex planes, as a comparison between local maxima 
/ max and local minima / mm of the order parameter shows: 



f(Vv 



/max f(y v ±ir/2epH, 



2V2 



exp 



c2ooy 



2a(?{T) 



< 1. 



(103) 



At lower temperatures, when 



H c 2oo (T) 



«C 2 (T) 
aCo 

vraC(T) 



«1, 



2 aft 



(104) 



(105) 



In this "2D" regime, the energy of the Josephson interlayer coupling is small relative to the 
intralayer condensation energy.0 The transition to the normal phase occurs owing mainly to 
pair-breaking by the intralayer super currents, and the order parameter is almost unperturbed 



by the vortex planes: 



f(y) « 1 



a 2 a( 2 (T) 
12 <o 
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cos (2epHy) 



(106) 



This expression should be compared with Eq. (|73|) for intermediate fields in the same 
temperature range ( |104| ). In the limit a = (no Josephson coupling), equation ( |105| ) goes 
over into the well familiar one for an isolated thin superconducting film.0@ll! Equation 
( |105| ) explains the origin of the well-knownlil unphysical "low-temperature" divergence of 
the LD model: taking a formal limit a — ► while keeping a( 2 (T) / a£ =const, we get 
H c2o o{T) -> oo. 

Aside from microscopically determined parameters, dependence (|102|) for layered super- 
conductors was first obtained within the framework of the LD model.iiUl! Expressions of 
the type ( 105| ) were derived phenomenologically in Refs. (§0). In all these publications 



relations (0), (|3lD were implicitly adopted as physically plausible assumptions. The very 
fact that these results are contained as limiting cases in our Eqs. (f42|)-(f48|) once again 
demonstrates the generality and self-consistency of the approach of this paper. 

Finally, we emphasize that the concept of Josephson vortex planes applies both in lim- 
its (|10l ) and (104). Contrary to previous suggestions,il@ there is no transition from the 



"Abrikosov-core regime" to the " Josephson-core regime" at T*\ The existence of Abrikosov 
vortices with normal cores in the thin-layer limit is not allowed by the solutions of Eq. (|9l|) . 
[Mathematically, the function ceo(t,q) is strictly positive.] 



F. Size-effects: Oscillations of T c w 

Aside from a special case epHL = 7ck/2 ($ = k& , k = 0, 1, 2, . . .), for multilayers with 
finite S-layer length W = 2L only approximate solutions of the boundary problem (|9T|). 
(P?2]) can be obtained. Using Galerkin's method,il we have found two groups of solutions 
corresponding to the smallest eigenvalues [A] c : 
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fn,N v =2n+i(t) oc cosh (/it) ce (t, q), fn,N v =2n(t) oc cosh (/it) ce ( - - t, q ) , (107) 



7T 



/i = coth (fiepHL) 



^(epHL,q) 



ce (epHL, q) 



(108) 



[A(T,if)] c = a (g)-P 



(109) 



and 



/«/,jv„=2n+i(*) oc cos (ut) ce (t, g), /„,jv„=2n(*) oc cos (i/i) ce - t, 



(110) 



i/ = — cot (uepHL) 



^(e P HL,q) 



ce (epHL, q) 



(111) 



[A(T,if)] c = a (g) + z/ 



112) 



where Eqs. ( |108| ), ( |111| ) implicitly define parameters /i and i/, and Eqs. ( |109| ), (|112|) deter- 
mine the critical point. 

From Eqs. ( |109| ), ( |112| ) we see, that the eigenvalue corresponding to the eigenfunctions / M 
is smaller than ao(q), while that corresponding to the eigenfunctions f v is larger. Physically, 
this means that with / M in a finite multilayer we can achieve higher values of the critical 
parameters T c and H c2 than in an infinite one. [Compare with Eq. (|93|).1 At epHL = irk/2 
($ = k&o, k = 0, 1, 2, . . .), these equations yield /i = v — 0. 

For epHL — > oo, /i is a bounded, oscillating function of epHL and does not tend to any 
limit. On the contrary, z/ — > 7c/2epHL, when epHL —>■ oo. This signifies that at certain 
values of epHL the solution / M becomes unstable and gives way to the solution /„ with lower 
values of the critical parameters, presumably by means of a first-order phase transition. As 
the parameters /i and v can enter the free energy only via the combinations fiepHL and 
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uepHL, we expect the transitions / M <-> to occur when \p,epHL]^ = [uepHL]^ , hence the 
relation 

cot \p,epHL]^ = — coth [p,epHL]^ 



with the numerical solution [/lepHL]^ ~ 2.37. Thus, the solution with ( |109| ) is realized 
when 



epHL ( ep #L, g) 



while for 



ce (epHL, q) 



epHL 



< [pepHL]janh [piepHL]^ w 2.33, 



(113) 



^(epHL^q) 



ce (epHL, q) 



> 2.33, 



(114) 



the system "chooses" /„ with (112). The condition ( |113[ ) is met, for instance, when <I> = 
pWH ps fc$ = 2L, A; = 0, 1, 2, . . .). Because of the oscillating character of the left-hand 
sides of Eqs. ( |113| ), ( 114] ), the system oscillates between the states with f ^ and /„ with 
increasing epHL. For larger epHL, the domain of existence of / M becomes narrower, while 
that of f v widens. For epHL — > oo, the solution /„ goes over smoothly into that of an 
infinite multilayer. 

We want to point out here that the exact character of the transitions «-> f v can only 
be established by solving the nonlinear boundary problem and comparing the corresponding 
free energies, which is beyond the scope of the present paper. 

As an important application of the above results, we consider the critical temperature 
of a finite multilayer T c \y in the field range given by (^H), and with W <C p\j/X(T). Under 
such circumstances, the condition ( |113j ) is satisfied, and the solution of Eq. ( |108| ) is 



a 



o 



a£ (epH) vr$ 
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which on substituting into Eq. ( |109| ) yields: 
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(115) 



Thus, in a finite-size multilayer the critical temperature can exhibit the same oscillations 
with changing the flux through the "elementary cell" as the total Josephson current / does 
[see Eq. fl80|) 1. However, a significant difference lies in the fact that while the oscillations of 
I are observable in any types of Josephson systems, the oscillations of T c is a novel feature, 
because in a single Josephson junction with thick superconducting electrodes any shifts of 
T c are negligible. In the limit $ 3> $o, equation ( |115| ) reduces to Eq. (|^), as anticipated. 



IV. DISCUSSION 

Based solely on the microscopic Hamiltonian (JI|), we have constructed a self-consistent 
theory that provides a comprehensive, unified picture of physical effects in S/I multilayers 
in parallel magnetic fields in the GL regime. 

Employing rigorous technique of variational calculus, we have derived in Sec. II funda- 
mental constraint relations (|27|), (|28D and solved a nontrivial problem of exact minimization 
of the microscopic functional (||). Up until the present study, such a problem has not been 
solved even for the much simpler phenomenological LD functional. Surprisingly, even mu- 
tual dependence and incompleteness of the Euler-Lagrange equations for <p n and A were not 
noticed in previous publications. This incompleteness fully manifests itself in unphysical 
degrees of freedom and an irrelevant length scale of equations for the phase differences pro- 
posed, e. g., in Refs. (®@): Making use of constraints of the type (|2lf), one can reduce these 
equations to our Eq. ( [415] ) with f(y) = 1. Emerging as a direct mathematical consequence 
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of such general physical properties as gauge invariance and Josephson interlayer coupling, 
constraints (0), (p8|) should apply to any superconducting weakly coupled periodic struc- 
ture. The discovery of their role in minimizing the free energy makes further progress in the 
development of the theory possible. 

In the thin-layer limit which corresponds to the domain of validity of the phenomenologi- 
cal LD model, we have derived a remarkably simple, closed set of self-consistent microscopic 
mean-field equations (f42D-(fi~8|) and the generating functional (f4~9D . The fact that the solution 



(|67|), (|70|) and (|73|) of these equations describing the vortex state in an infinite multilayer 
reproduces the result obtained by Theodorakis0 in the framework of the LD model is not 
an occasional coincidence. The application of the mathematical methods of this paper al- 
lows us to obtain the complete exact solution of the LD model in parallel fields as well. 
(This solution will be published elsewhere.) The resulting mean-field equations are merely 
a limiting case of our Eqs. (P^)-(]4"5|) for a/p — > 0. As our equations contain more physical 
information, we propose that they should replace the LD model in parallel fields. 

Concerning some major physical results of Sec. Ill in the thin-layer limit, we have 
completely revised previous calculations!!! of H c i based on an invalid assumption of single 
Josephson vortex penetration and refuted the concept! of a triangular Josephson vortex 
lattice. Our consideration envisages simultaneous and coherent penetration in the form of 
the "vortex planes" . Our prediction of the superheating field H s for semi-infinite multilayers 
implies hysteretic behavior of the magnetization. In the vortex state, the magnetization 
should exhibit jumps due to the vortex-plane penetration. We have fully clarified the wide- 
spread!S misunderstanding of the physics of the Fraunhofer oscillations: our self-consistent 
treatment of the Josephson effect unambiguously proves that the Fraunhofer pattern occurs 
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due to successive penetration of the vortex planes and their pinning by the edges of the 
sample. Our prediction of novel size-effects in finite multilayers, a series of first-order phase 
transitions to the normal state and oscillations of the critical temperature versus the applied 
field, should stimulate further experimental investigation. 

The results of our investigation directly apply to artificial superconductor/insulatori! 
and superconductor/semiconductoril'ii'i multilayers. As regards the high-T c superconduc- 
tors BSCCO and TBCCO, believed to be atomic-scale weakly coupled superlattices,0 the 
application is restricted by the limitation (f|). However, we expect that such basic features 
of the thin-layer limit as simultaneous and coherent penetration in the form of the vortex 
planes will hold. For high-T c samples exhibiting a clear Fraunhofer pattern,!!] we anticipate 
the presence of the related effect of oscillations of the critical temperature ( |115|) as well. 



As to direct experimental verification of basic concepts of our theory, the best evidence 
is provided by the recent magnetization and polarized neutron reflectivity measurements 
on Nb/Si multilayers in parallel fields.^ These measurements clearly revealed simultaneous 
penetration of Josephson vortices into all Si layers and a companion effect of jumps of the 
magnetization, exactly as predicted in our paper. The distribution of the magnetic field 
attributed to Josephson vortices was found to have the periodicity of the Nb/Si layering, 
in agreement with the general consideration of Sec. II. A. (The experimental conditions^ 
did not fully match the requirements of the thin-layer limit for which the screening by the 
intralayer currents could be neglected.) Finally, it is quite natural that our general implicit 
expression ( |10U| ) for H c2oo (T) exhibits the so-called "3D-2D crossover", well-known from the 
experiment,! and is free from the unphysical "low-temperature" divergence, typicali of the 
phenomenological LD model. 
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APPENDIX A: THE APPLICATION OF MATHIEU FUNCTIONS 



The canonical form of the Mathieu equation H 1 

d 2 f 



dt 2 

If f(t) is a solution to (|A1|) , then / f~ — t) satisfies 



+ (a - 2gcos2t) / = 0. (Al) 



+ (a + 2qcos2t)f = 0. (A2) 

In the class of periodic solutions of (|A"T|), the smallest eigenvalue a^q) is a non-positive, 
continuous, even, monotonously decreasing function of q. The corresponding eigenfunction 
ceo(£, q) has a period n, is even and strictly positive. 

For <q <C 1, we have the asymptotics 



ce (t,g) w 



1 - | cos It + . . 



(A3) 



a (q) « -£ + . . . (A4) 



For g 3> 1, 



a (g)~2^-2g, (A5) 



but there is no uniform asymptotics for ceo(t,q). In this case, the behavior of ce (t,q) may 
be characterized by the formulas 
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Ce (t,g)~f-J gV8 e -V5cos 2 t/4 j |cos£| < 



2 l/4 

7r 



ce (0, g) ~ 2v / 2e- 2 v / ^ce (7r/2, g) ~ 2 (2tt) 1/4 g 1/8 e - 2 A 
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FIGURES 



Fig. 1. Geometry of the problem. Alternating superconducting layers and nonsupercon- 
ducting barriers are shown by white and grey rectangles, respectively. The system is 
supposed to be infinite in the x and z directions. An external magnetic field H is 
applied along the z axis. 

Fig. 2. Vortex state in a finite LD multilayer (in cross section). Josephson vortices (i. 
e., singularities of the amplitude of condensation) are conventionally denoted by black 
dots. The "vortex planes" [i. e., maxima of the microscopic magnetic field h(y)] are 
shown by dashed lines. Arrows show the direction of supercurrents. 
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